deficit_Model_Plots

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(knitr)
library(ggplot2)
library(dplyr)
library(lubridate)
library(MMWRweek)
library(gridExtra)

## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine

library(readxl)
library(excessmort)
library(plotly)

## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout

library(patchwork)
library(scales)

## 
## Attaching package: 'scales'
## 
## The following object is masked from 'package:purrr':
## 
##     discard
## 
## The following object is masked from 'package:readr':
## 
##     col_factor

library(cowplot)

## 
## Attaching package: 'cowplot'
## 
## The following object is masked from 'package:patchwork':
## 
##     align_plots
## 
## The following object is masked from 'package:lubridate':
## 
##     stamp

library(broom)
library(gt)

## Warning: Can't find generic `as.gtable` in package gtable to register S3 method.
## ℹ This message is only shown to developers using devtools.
## ℹ Do you need to update gtable to the latest version?

## 
## Attaching package: 'gt'
## 
## The following object is masked from 'package:cowplot':
## 
##     as_gtable

library(webshot2)
library(stringr)
library(purrr)
library(kableExtra)

## 
## Attaching package: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows

source("~/Desktop/Thesis Work/Hospitalization/11-26 data/Deficit_Model_1.0.R")

allcause_state_byweek <- read_excel("allcause_state_byweek_suppressed_25mar21.xlsx")
icd10categories <- read_excel("icd10categories_updated_321.xlsx")

#weekly state data
all_cause_state_byweek <- allcause_state_byweek %>%
  rename(weekly_count_statelevel = count)

all_cause_state_by_week <- all_cause_state_byweek %>%
  separate(mmwr_week, into = c("MMWRyear", "MMWRweek"), sep = "-", convert = TRUE) %>%
  mutate(
    date = MMWRweek2Date(MMWRyear = MMWRyear, MMWRweek = MMWRweek, MMWRday = 1),
    statecount_numeric = as.numeric(ifelse(weekly_count_statelevel == "<5", NA, weekly_count_statelevel)),
    statecount_suppressed = weekly_count_statelevel == "<5",
    category = as.character(category)
  )

all_cause_state_by_week %>%
  group_by(category) %>%
  reframe(missing_total = sum(statecount_suppressed), 
            count_total = sum(statecount_numeric, na.rm = TRUE))

## # A tibble: 10 × 3
##    category               missing_total count_total
##    <chr>                          <int>       <dbl>
##  1 cardiovascular                     0      373147
##  2 covid                              1       58262
##  3 digestive                          0      182701
##  4 infectious                         0      112279
##  5 infectious respiratory             0       29816
##  6 injury                             0      331205
##  7 other                              0     1740863
##  8 pneumonia                          0       97954
##  9 renal                              0      102199
## 10 respiratory                        0      164424

all_cause_state_by_week

## # A tibble: 3,140 × 7
##    MMWRyear MMWRweek category               weekly_count_statelevel date      
##       <int>    <int> <chr>                  <chr>                   <date>    
##  1     2019        1 cardiovascular         843                     2018-12-30
##  2     2019        1 digestive              441                     2018-12-30
##  3     2019        1 infectious             299                     2018-12-30
##  4     2019        1 infectious respiratory 240                     2018-12-30
##  5     2019        1 injury                 691                     2018-12-30
##  6     2019        1 other                  4090                    2018-12-30
##  7     2019        1 pneumonia              427                     2018-12-30
##  8     2019        1 renal                  220                     2018-12-30
##  9     2019        1 respiratory            522                     2018-12-30
## 10     2019        1 covid                  0                       2018-12-30
## # ℹ 3,130 more rows
## # ℹ 2 more variables: statecount_numeric <dbl>, statecount_suppressed <lgl>

str(icd10categories)

## tibble [869 × 5] (S3: tbl_df/tbl/data.frame)
##  $ icd10code          : chr [1:869] "A00" "A01" "A02" "A03" ...
##  $ category_original  : chr [1:869] "infectious" "infectious" "infectious" "infectious" ...
##  $ category_1113update: chr [1:869] "infectious" "infectious" "infectious" "infectious" ...
##  $ category_123update : chr [1:869] "infectious" "infectious" "infectious" "infectious" ...
##  $ category_0321update: chr [1:869] "infectious" "infectious" "infectious" "infectious" ...

names(icd10categories)

## [1] "icd10code"           "category_original"   "category_1113update"
## [4] "category_123update"  "category_0321update"

dim(icd10categories)

## [1] 869   5

summary(icd10categories)

##   icd10code         category_original  category_1113update category_123update
##  Length:869         Length:869         Length:869          Length:869        
##  Class :character   Class :character   Class :character    Class :character  
##  Mode  :character   Mode  :character   Mode  :character    Mode  :character  
##  category_0321update
##  Length:869         
##  Class :character   
##  Mode  :character

sapply(icd10categories, class)

##           icd10code   category_original category_1113update  category_123update 
##         "character"         "character"         "character"         "character" 
## category_0321update 
##         "character"

glimpse(icd10categories)

## Rows: 869
## Columns: 5
## $ icd10code           <chr> "A00", "A01", "A02", "A03", "A04", "A05", "A06", "…
## $ category_original   <chr> "infectious", "infectious", "infectious", "infecti…
## $ category_1113update <chr> "infectious", "infectious", "infectious", "infecti…
## $ category_123update  <chr> "infectious", "infectious", "infectious", "infecti…
## $ category_0321update <chr> "infectious", "infectious", "infectious", "infecti…

colSums(is.na(icd10categories))

##           icd10code   category_original category_1113update  category_123update 
##                   0                 418                 418                   1 
## category_0321update 
##                   0

str(allcause_state_byweek)

## tibble [3,140 × 3] (S3: tbl_df/tbl/data.frame)
##  $ mmwr_week: chr [1:3140] "2019-01" "2019-01" "2019-01" "2019-01" ...
##  $ category : chr [1:3140] "cardiovascular" "digestive" "infectious" "infectious respiratory" ...
##  $ count    : chr [1:3140] "843" "441" "299" "240" ...

names(allcause_state_byweek)

## [1] "mmwr_week" "category"  "count"

dim(allcause_state_byweek)

## [1] 3140    3

summary(allcause_state_byweek)

##   mmwr_week           category            count          
##  Length:3140        Length:3140        Length:3140       
##  Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character

sapply(allcause_state_byweek,class)

##   mmwr_week    category       count 
## "character" "character" "character"

glimpse(allcause_state_byweek)

## Rows: 3,140
## Columns: 3
## $ mmwr_week <chr> "2019-01", "2019-01", "2019-01", "2019-01", "2019-01", "2019…
## $ category  <chr> "cardiovascular", "digestive", "infectious", "infectious res…
## $ count     <chr> "843", "441", "299", "240", "691", "4090", "427", "220", "52…

colSums(is.na(allcause_state_byweek))

## mmwr_week  category     count 
##         0         0         0

Deficit Function Used

deficit_plot <- function(fit, 
                        title     = "", 
                        ylim      = NULL,
                        show.data = TRUE,
                        alpha     = 0.05){

  requireNamespace("ggplot2")

  z <- qnorm(1 - alpha/2)


  p <- with(fit, data.frame(date = date,
                        y        = (observed - expected)/expected,
                        increase = fitted,
                        sd       = sd,
                        se       = se)) %>%
    ggplot(aes(date, y)) +
    geom_ribbon(aes(ymin = increase - z * se, ymax = increase + z * se), alpha = 0.5, fill = "#3366FF") + #needs to be edited
    geom_hline(yintercept = 0) +
    ylab("% Change from expected count") +
    scale_y_continuous(labels = scales::percent) +
    scale_x_date(date_labels = "%b %Y") +
    ggtitle(title) +
    xlab("Date") 

  if(show.data) p <- p + geom_point(alpha = 0.3)

  if(!is.null(ylim)) p <- p + coord_cartesian(ylim = ylim)

  return(p + geom_line(aes(y = increase), size=0.70, col="#3366FF"))
}

deficit_cumulative <- function(fit, start, end){
  if (!"curve_fit" %in% attr(fit, "type"))
    stop("This is not the correct deficit_model fit, needs curve fit.")

  ind             <- which(fit$date %in% seq(start, end, by = "day"))
  n               <- length(ind) # returns the # of elements in the ind object
  A               <- matrix(1, n, n)
  A[upper.tri(A)] <- 0                                         # It sets all the upper triangular elements of matrix A (excluding the diagonal) to 0
  A               <- sweep(A, 2, fit$expected[ind], FUN = "*") # This line multiplies each column (2 means columns) of the matrix A by a corresponding value from fit$expected[ind]
                                                               # So, now it's scaling each column of A by expected values from the model fit (fit$expected[ind])
  fhat <- matrix(fit$fitted[ind], ncol = 1)                    # creates a column vector (matrix with 1 column) called fhat from the values in fit$fitted[ind]

 
   # Below, this code is constructing a data frame with:
  
  fit_deficit   <- A %*% fhat                                    # multiply design matrix A by fitted values fhat to get cumulative modeled effect.
  obs_deficit   <- cumsum(fit$observed[ind] - fit$expected[ind])
  varcov       <- fit$x[ind,] %*% fit$betacov %*% t(fit$x[ind,]) # Compute variance-covariance matrix of the linear predictor based on model design matrix x and estimated betacov
  diag(varcov) <- fit$se[ind]^2                                  # Replace diagonal entries of varcov with squared standard errors — probably to correct for over/under-dispersion or model misspecification.
  fitted_se    <- sqrt(diag(A %*%  varcov %*% t(A)))
  sd           <- sqrt(diag(A %*% (fit$cov[ind,ind] + diag(fit$log_expected_se[ind]^2)) %*% t(A))) # this computes the variance of the observed cumulative deficit.
                                                                                                   # It adds fit$cov (covariance of observed data) and extra uncertainty from fit$log_expected_se.
  data.frame(date     = fit$date[ind], # The time points
             observed = obs_deficit,   # Cumulative observed - expected difference
             sd       = sd,            # Modeled cumulative deficit/excess from the model.
             fitted   = fit_deficit,   # Standard deviation of the observed cumulative deficit.
             se       = fitted_se)     # Standard error of the fitted cumulative deficit.
}

Deficit Model fits for the different categories

Cardio

exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 2, 20), 
                     by = "day")

# cardiovascular data
cardio <- all_cause_state_by_week %>%
  filter(category == "cardiovascular") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)
3

## [1] 3

# deficit model with COVID-19 dates excluded
cardio_deficit <- compute_expected(cardio, 
                               exclude = exclude_dates,
                               frequency = 52)

## Overall death rate is 8.85.

cardio_deficit

## # A tibble: 312 × 12
##    MMWRyear MMWRweek category       weekly_count_statelevel date      
##       <int>    <int> <chr>          <chr>                   <date>    
##  1     2019        2 cardiovascular 1233                    2019-01-06
##  2     2019        3 cardiovascular 1160                    2019-01-13
##  3     2019        4 cardiovascular 1201                    2019-01-20
##  4     2019        5 cardiovascular 1157                    2019-01-27
##  5     2019        6 cardiovascular 1286                    2019-02-03
##  6     2019        7 cardiovascular 1261                    2019-02-10
##  7     2019        8 cardiovascular 1141                    2019-02-17
##  8     2019        9 cardiovascular 1156                    2019-02-24
##  9     2019       10 cardiovascular 1208                    2019-03-03
## 10     2019       11 cardiovascular 1232                    2019-03-10
## # ℹ 302 more rows
## # ℹ 7 more variables: statecount_numeric <dbl>, statecount_suppressed <lgl>,
## #   population <dbl>, outcome <dbl>, log_expected_se <dbl>, expected <dbl>,
## #   excluded <lgl>

cardio_deficit <- cardio_deficit %>%
  mutate(date = as.Date(paste(MMWRyear, MMWRweek, 1, sep = "-"), "%Y-%U-%u"))

## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `date = as.Date(paste(MMWRyear, MMWRweek, 1, sep = "-"),
##   "%Y-%U-%u")`.
## Caused by warning in `strptime()`:
## ! (0-based) yday 369 in year 2020 is invalid

combined_plot_1 <- ggplot(cardio_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Cardiovascular Hospitalizations in MA",
       x = "Date",
       y = "Weekly Hospitilizations",
       color = "Data Point Status") +
  theme_minimal()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

combined_plot_1

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("combined_plot_1.png", plot = combined_plot_1, width = 8, height = 6, dpi = 300)

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

fit_x <- cardio_deficit %>% 
  deficit_model(exclude = exclude_dates,
               start = as.Date("2019-01-14"),
               end = as.Date("2024-12-31"),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_x <- deficit_plot(fit_x, 
            title = "Deviations from Expected Counts for Cardio in MA"
)
deficit_plot_x

fit_x$detected_intervals

##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-09-23 2019-12-02       8.144139       8.669496    0.07650853    12062
## 2 2020-01-20 2020-02-03       8.276251       8.703193    0.14678711     3343
## 3 2020-03-02 2021-02-15       7.265876       8.874462    0.03594969    49893
## 4 2021-12-27 2022-01-31       8.304721       9.024090    0.10569035     6709
## 5 2023-05-01 2023-05-01       9.677494       9.573722    0.26665497     1303
##    expected     deficit        sd
## 1 12840.088   -778.0880 113.31411
## 2  3515.453   -172.4533  59.29126
## 3 60938.774 -11045.7736 246.85780
## 4  7290.145   -581.1446  85.38234
## 5  1289.028     13.9721  35.90303

# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_1 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-Cardio",
       title    = "Cumulative Deficit Hospitilization for Cardio Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_1

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_1_adj <- combined_plot_1 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_x_adj <- deficit_plot_x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

cumulative_deficit_plot_1_adj <- cumulative_deficit_plot_1 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_1 <- combined_plot_1_adj /
                       deficit_plot_x_adj /
                       cumulative_deficit_plot_1_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_1

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("Cardio_combined_figure.png",
       plot = final_combined_plot_1,
       width = 11,      
       height = 10.5,    
       dpi = 300,
       bg = "white")

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("cardio_combined_figure.png",
       plot = final_combined_plot_1,
       width = 11,       # good for standard document width
       height = 10.5,    # leaves space for three stacked plots
       dpi = 300,
       bg = "white")

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

Digestive

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 2, 27), 
                     by = "day")

# digestive data
digestive <- all_cause_state_by_week %>%
  filter(category == "digestive") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
digestive_deficit <- compute_expected(digestive, 
                               exclude = exclude_dates,
                               frequency = 52)

## Overall death rate is 4.33.

digestive_deficit

## # A tibble: 312 × 12
##    MMWRyear MMWRweek category  weekly_count_statelevel date      
##       <int>    <int> <chr>     <chr>                   <date>    
##  1     2019        2 digestive 535                     2019-01-06
##  2     2019        3 digestive 570                     2019-01-13
##  3     2019        4 digestive 553                     2019-01-20
##  4     2019        5 digestive 576                     2019-01-27
##  5     2019        6 digestive 609                     2019-02-03
##  6     2019        7 digestive 599                     2019-02-10
##  7     2019        8 digestive 602                     2019-02-17
##  8     2019        9 digestive 562                     2019-02-24
##  9     2019       10 digestive 550                     2019-03-03
## 10     2019       11 digestive 622                     2019-03-10
## # ℹ 302 more rows
## # ℹ 7 more variables: statecount_numeric <dbl>, statecount_suppressed <lgl>,
## #   population <dbl>, outcome <dbl>, log_expected_se <dbl>, expected <dbl>,
## #   excluded <lgl>

# Plotting
combined_plot_2 <- ggplot(digestive_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly digestive ED Hospitilizations in MA",
       x = "Date",
       y = "Weekly Hospitilizations",
       color = "Data Point Status") +
  theme_minimal()

combined_plot_2

fit_2x <- digestive_deficit %>% 
  deficit_model(exclude = exclude_dates,
               start = min(.$date),
               end = max(.$date),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_2x <- deficit_plot(fit_2x, 
            title = "Deviations from Expected Counts for Digestive in MA"
)
deficit_plot_2x

fit_2x$detected_intervals

##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-09-15 2019-10-13       3.986860       4.261767    0.07956446     2684
## 2 2020-03-01 2020-09-06       3.519113       4.442692    0.03432838    13267
## 3 2020-10-11 2021-03-14       3.714835       4.193817    0.03680019    11504
## 4 2021-12-05 2022-01-30       3.759757       4.186787    0.05877984     4556
## 5 2022-05-01 2022-05-22       4.368984       4.646495    0.09288423     2353
## 6 2022-09-18 2022-10-02       4.263148       4.543836    0.10606203     1722
##    expected    deficit        sd
## 1  2869.071  -185.0706  53.56371
## 2 16748.877 -3481.8769 129.41745
## 3 12987.299 -1483.2986 113.96183
## 4  5073.467  -517.4671  71.22827
## 5  2502.459  -149.4588  50.02458
## 6  1835.377  -113.3774  42.84131

# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_2x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_2 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-Digestive",
       title    = "Cumulative Deficit Hospitilization for Digestive Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_2

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_2_adj <- combined_plot_2 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_2x_adj <- deficit_plot_2x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

cumulative_deficit_plot_2_adj <- cumulative_deficit_plot_2 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_2 <- combined_plot_2_adj /
                       deficit_plot_2x_adj /
                       cumulative_deficit_plot_2_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_2

ggsave("Digestive_combined_figure.png",
       plot = final_combined_plot_2,
       width = 11,      
       height = 10.5,    
       dpi = 300,
       bg = "white")

Infectious

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 6, 5), 
                     by = "day")

# infectious data
infectious <- all_cause_state_by_week %>%
  filter(category == "infectious") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# deficit model with COVID-19 dates excluded
infectious_deficit <- compute_expected(infectious,
                               exclude = exclude_dates,
                               frequency = 52)

## Overall death rate is 2.66.

# Plotting
combined_plot_3 <- ggplot(infectious_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Infectious ED Hospitilizations in MA",
       x = "Date",
       y = "Weekly Hospitilizations",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_3

fit_3x <- infectious_deficit %>% 
  deficit_model(exclude = exclude_dates,
               start = min(.$date),
               end = max(.$date),
               knots.per.year = 12,
               verbose = FALSE)

deficit_plot_3x <- deficit_plot(fit_3x, 
            title = "Deviations from Expected Counts for Infectious in MA"
)
deficit_plot_3x

fit_3x$detected_intervals

##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2020-05-03 2020-08-02       2.384625       2.718529    0.03797626     4495
## 2 2020-10-04 2021-05-30       2.231734       2.726452    0.02405327    10517
## 3 2022-01-23 2022-02-27       2.530161       2.785967    0.05872481     2044
## 4 2022-12-25 2023-01-22       2.565316       2.820316    0.06472516     1727
## 5 2023-06-11 2023-07-02       2.519639       2.766860    0.07167585     1357
## 6 2024-10-20 2024-12-15       2.525210       2.768028    0.04779399     3060
##    expected    deficit        sd
## 1  5124.405  -629.4047  71.58495
## 2 12848.349 -2331.3491 113.35056
## 3  2250.654  -206.6542  47.44106
## 4  1898.669  -171.6689  43.57372
## 5  1490.146  -133.1456  38.60240
## 6  3354.243  -294.2430  57.91583

# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_3x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_3 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-Infectious",
       title    = "Cumulative Deficit Hospitilization for Infectious Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_3

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_3_adj <- combined_plot_3 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_3x_adj <- deficit_plot_3x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

cumulative_deficit_plot_3_adj <- cumulative_deficit_plot_3 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_3 <- combined_plot_3_adj /
                       deficit_plot_3x_adj /
                       cumulative_deficit_plot_3_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_3

ggsave("Infectious_combined_figure.png",
       plot = final_combined_plot_3,
       width = 11,       
       height = 10.5,    
       dpi = 300,
       bg = "white")

infectious_respiratory

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 6, 12), 
                     by = "day")

# infectious respiratory data
infectious_respiratory <- all_cause_state_by_week %>%
  filter(category == "infectious respiratory") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
infectious_respiratory_deficit <- compute_expected(infectious_respiratory,
                               exclude = exclude_dates,
                               frequency = 52)

## Overall death rate is 0.699.

# Plotting
combined_plot_4 <- ggplot(infectious_respiratory_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Infectious Respiratory ED Hospitilizations in MA",
       x = "Date",
       y = "Weekly Hospitilizations",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_4

# -- Fitting excess model
fit_4x <- infectious_respiratory_deficit %>% 
  deficit_model(exclude = exclude_dates,
                start = min(.$date),
               end = max(.$date),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_4x <- deficit_plot(fit_4x, 
            title = "Deviations from Expected Counts for Infectious Respiratory in MA"
)
deficit_plot_4x

# -- Intervals of inordinate mortality found by the excess model
fit_4x$detected_intervals

##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-10-13 2019-12-01      0.6340876      1.2580696    0.03417564      683
## 2 2020-03-29 2021-05-30      0.1241042      0.8568669    0.01013141     1036
## 3 2021-10-31 2022-03-20      0.4714432      1.3271874    0.02166535     1333
## 4 2023-01-15 2023-02-26      0.6535837      1.2568241    0.03651720      616
## 5 2023-04-09 2023-04-23      0.4208683      0.7014588    0.04167251      170
## 6 2024-11-03 2024-11-10      0.6424430      1.0037101    0.06105179      173
##    expected     deficit       sd
## 1 1355.1149  -672.11494 36.81189
## 2 7152.9723 -6116.97228 84.57525
## 3 3752.6065 -2419.60646 61.25852
## 4 1184.5517  -568.55175 34.41732
## 5  283.3380  -113.33805 16.83265
## 6  270.2837   -97.28366 16.44031

# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_4x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

cumulative_deficit_plot_4 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-infectious respiratory",
       title    = "Cumulative Deficit Hospitilization for infectious respiratory Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_4

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_4_adj <- combined_plot_4 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_4x_adj <- deficit_plot_4x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

cumulative_deficit_plot_4_adj <- cumulative_deficit_plot_4 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_4 <- combined_plot_4_adj /
                       deficit_plot_4x_adj /
                       cumulative_deficit_plot_4_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_4

ggsave("Infectious Respiratory_combined_figure.png",
       plot = final_combined_plot_4,
       width = 11,       # good for standard document width
       height = 10.5,    # leaves space for three stacked plots
       dpi = 300,
       bg = "white")

injury

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 4, 10), 
                     by = "day")

# injury data
injury <- all_cause_state_by_week %>%
  filter(category == "injury") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
injury_deficit <- compute_expected(injury,
                               exclude = exclude_dates,
                               frequency = 52)

## Overall death rate is 7.86.

# Plotting
combined_plot_5 <- ggplot(injury_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Injury ED Hospitilizations in MA",
       x = "Date",
       y = "Weekly Hospitilizations",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_5

# -- Fitting excess model
fit_5x <- injury_deficit %>%                                                     
  deficit_model(exclude = exclude_dates,
                start = min(.$date),
               end = max(.$date),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_5x <- excess_plot(fit_5x, 
            title = "Deviations from Expected Counts for Injuries in MA"
)
deficit_plot_5x

# -- Intervals of inordinate mortality found by the excess model
fit_5x$detected_intervals

##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-10-20 2019-11-24       6.756174       7.279676    0.09492698     5458
## 2 2020-01-12 2020-09-13       6.349953       7.623419    0.03965819    30779
## 3 2020-10-04 2021-02-07       6.578836       7.492868    0.05411984    16830
## 4 2021-03-07 2021-04-11       7.120101       7.580370    0.09686766     5752
## 5 2021-05-02 2021-05-16       7.600386       7.944070    0.14023944     3070
## 6 2021-07-25 2021-08-08       7.979167       8.372155    0.14396843     3223
## 7 2021-12-19 2022-01-16       7.140402       7.682093    0.10682282     4807
## 8 2022-07-31 2022-08-21       8.268205       8.668967    0.12687117     4453
##    expected    deficit        sd
## 1  5880.914  -422.9136  76.68712
## 2 36951.644 -6172.6438 192.22810
## 3 19168.280 -2338.2799 138.44956
## 4  6123.830  -371.8301  78.25490
## 5  3208.823  -138.8232  56.64648
## 6  3381.738  -158.7384  58.15272
## 7  5171.673  -364.6727  71.91434
## 8  4668.838  -215.8381  68.32890

# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_5x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_5 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-Injury",
       title    = "Cumulative Deficit Hospitilization for Injury Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_5

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_5_adj <- combined_plot_5 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_5x_adj <- deficit_plot_5x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

cumulative_deficit_plot_5_adj <- cumulative_deficit_plot_5 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_5 <- combined_plot_5_adj /
                       deficit_plot_5x_adj /
                       cumulative_deficit_plot_5_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_5

ggsave("Injury_combined_figure.png",
       plot = final_combined_plot_5,
       width = 11,       # good for standard document width
       height = 10.5,    # leaves space for three stacked plots
       dpi = 300,
       bg = "white")

pneumonia

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2022, 5, 7), 
                     by = "day")

# pneumonia data
pneumonia <- all_cause_state_by_week %>%
  filter(category == "pneumonia") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
pneumonia_deficit <- compute_expected(pneumonia,
                               exclude = exclude_dates,
                               frequency = 52)

## Overall death rate is 2.32.

# Plotting
combined_plot_6 <- ggplot(pneumonia_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Pneumonia ED Hospitilizations in MA",
       x = "Date",
       y = "Weekly Hospitilizations",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_6

# -- Fitting excess model
fit_6x <- pneumonia_deficit %>% 
  deficit_model(exclude = exclude_dates,
               start = min(.$date),
               end = max(.$date),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_6x <- deficit_plot(fit_6x, 
            title = "Deviations from Expected Counts for Pneumonia in MA"
)
deficit_plot_6x

# -- Intervals of inordinate mortality found by the excess model
fit_6x$detected_intervals

##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-10-27 2019-12-01       2.532637       2.967751    0.06061043     2046
## 2 2020-04-26 2021-07-04       1.489661       2.618500    0.01756974    12636
## 3 2021-09-05 2022-05-22       1.947460       2.813849    0.02345135     9964
## 4 2022-07-17 2022-08-28       1.694437       2.049678    0.04663405     1597
## 5 2022-11-06 2022-11-13       2.536350       2.922000    0.10416801      683
## 6 2023-01-22 2023-02-19       2.407862       2.978610    0.06651675     1621
## 7 2023-04-23 2023-05-21       2.294970       2.662542    0.06288869     1545
## 8 2024-01-28 2024-02-11       2.619286       2.974941    0.08581986     1058
##     expected    deficit        sd
## 1  2397.5084  -351.5084  48.96436
## 2 22211.3287 -9575.3287 149.03466
## 3 14396.7975 -4432.7975 119.98666
## 4  1931.8138  -334.8138  43.95240
## 5   786.8496  -103.8496  28.05084
## 6  2005.2341  -384.2341  44.77984
## 7  1792.4540  -247.4540  42.33738
## 8  1201.6586  -143.6586  34.66495

# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_6x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_6 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-pneumonia",
       title    = "Cumulative Deficit Hospitilization for pneumonia Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_6

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_6_adj <- combined_plot_6 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_6x_adj <- deficit_plot_6x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

cumulative_deficit_plot_6_adj <- cumulative_deficit_plot_6 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_6 <- combined_plot_6_adj /
                       deficit_plot_6x_adj /
                       cumulative_deficit_plot_6_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_6

ggsave("Pneumonia_combined_figure.png",
       plot = final_combined_plot_6,
       width = 11,       
       height = 10.5,    
       dpi = 300,
       bg = "white")

Renal

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 2, 27), 
                     by = "day")

# renal data
renal <- all_cause_state_by_week %>%
  filter(category == "renal") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
renal_deficit <- compute_expected(renal,
                               exclude = exclude_dates,
                               frequency = 52)

## Overall death rate is 2.42.

# Plotting
combined_plot_7 <- ggplot(renal_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Renal ED Hospitilizations in MA",
       x = "Date",
       y = "Weekly Hospitilizations",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_7

# -- Fitting excess model
fit_7x <- renal_deficit %>% 
  deficit_model(exclude = exclude_dates,
                start = min(.$date),
               end = max(.$date),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_7x <- deficit_plot(fit_7x, 
            title = "Deviations from Expected Counts for Renal in MA"
)
deficit_plot_7x

# -- Intervals of inordinate mortality found by the excess model
fit_7x$detected_intervals

##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-11-10 2019-12-01       1.914332       2.068933    0.06198012     1031
## 2 2020-03-01 2020-07-12       1.798840       2.309561    0.02928593     4844
## 3 2020-09-06 2020-10-04       2.170195       2.326802    0.05879006     1461
## 4 2020-11-01 2021-02-28       1.966115       2.218355    0.03025440     4765
## 5 2021-12-26 2022-01-02       2.213272       2.330685    0.09303280      596
## 6 2022-05-15 2022-05-29       2.364289       2.607537    0.08034593      955
##    expected     deficit       sd
## 1 1114.2633   -83.26328 33.38058
## 2 6219.2915 -1375.29147 78.86248
## 3 1566.4295  -105.42948 39.57814
## 4 5376.3201  -611.32014 73.32339
## 5  627.6176   -31.61763 25.05230
## 6 1053.2544   -98.25443 32.45388

# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_7x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_7 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-renal",
       title    = "Cumulative Deficit Hospitilization for renal Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_7

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_7_adj <- combined_plot_7 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_7x_adj <- deficit_plot_7x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

cumulative_deficit_plot_7_adj <- cumulative_deficit_plot_7 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_7 <- combined_plot_7_adj /
                       deficit_plot_7x_adj /
                       cumulative_deficit_plot_7_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_7

ggsave("Renal_combined_figure.png",
       plot = final_combined_plot_7,
       width = 11,       # good for standard document width
       height = 10.5,    # leaves space for three stacked plots
       dpi = 300,
       bg = "white")

Respiratory

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 6, 12), 
                     by = "day")

# respiratory data
respiratory <- all_cause_state_by_week %>%
  filter(category == "respiratory") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
respiratory_deficit <- compute_expected(respiratory,
                               exclude = exclude_dates,
                               frequency = 52)

## Overall death rate is 3.89.

# Plotting
combined_plot_8 <- ggplot(respiratory_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Respiratory ED Hospitilizations in MA",
       x = "Date",
       y = "Weekly Hospitilizations",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_8

# -- Fitting excess model
fit_8x <- respiratory_deficit %>% 
  deficit_model(exclude = exclude_dates,
               start = min(.$date),
               end = max(.$date),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_8x <- deficit_plot(fit_8x, 
            title = "Deviations from Expected Counts for Respiratory in MA"
)
deficit_plot_8x

# -- Intervals of inordinate mortality found by the excess model
fit_8x$detected_intervals

##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-07-21 2019-07-21       2.978262       3.366797    0.15813126      401
## 2 2019-10-27 2019-11-24       3.920017       4.324783    0.08015053     2639
## 3 2020-03-08 2021-05-30       2.675694       4.168126    0.02182343    23417
## 4 2021-11-21 2022-03-27       3.632627       4.441859    0.04166917     9293
## 5 2023-04-23 2023-05-07       4.186401       4.696742    0.10783183     1691
##     expected      deficit        sd
## 1   453.3133    -52.31326  21.29115
## 2  2911.4933   -272.49331  53.95826
## 3 36478.3943 -13061.39427 190.99318
## 4 11363.1786  -2070.17861 106.59821
## 5  1897.1403   -206.14030  43.55617

# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_8x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_8 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-respiratory",
       title    = "Cumulative Deficit Hospitilization for respiratory Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_8

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_8_adj <- combined_plot_8 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_8x_adj <- deficit_plot_8x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

cumulative_deficit_plot_8_adj <- cumulative_deficit_plot_8 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_8 <- combined_plot_8_adj /
                       deficit_plot_8x_adj /
                       cumulative_deficit_plot_8_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_8

ggsave("Respiratory_combined_figure.png",
       plot = final_combined_plot_8,
       width = 11,      
       height = 10.5,   
       dpi = 300,
       bg = "white")

other

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 6, 5), 
                     by = "day")

# other data
other <- all_cause_state_by_week %>%
  filter(category == "other") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
other_deficit <- compute_expected(other,
                               exclude = exclude_dates,
                               frequency = 52)

## Overall death rate is 41.3.

# Plotting
combined_plot_9 <- ggplot(other_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Other ED Hospitilizations in MA",
       x = "Date",
       y = "Weekly Hospitilizations",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_9

# -- Fitting excess model
fit_9x <-  other_deficit %>% 
  deficit_model(exclude = exclude_dates,
                start = min(.$date),
               end = max(.$date),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_9x <- deficit_plot(fit_9x, 
            title = "Deviations from Expected Counts for Other in MA"
)
deficit_plot_9x

# -- Intervals of inordinate mortality found by the excess model
fit_9x$detected_intervals

##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2020-03-01 2021-05-23       35.68635       41.60681    0.06895008   312318
## 2 2021-08-22 2021-10-10       40.63174       42.79761    0.19933061    43766
## 3 2021-11-21 2022-02-13       36.00309       40.26895    0.15167805    63018
## 4 2022-03-06 2022-04-03       41.10001       42.97170    0.25264781    27669
## 5 2022-04-24 2022-05-08       42.59682       44.33594    0.33130392    17206
##    expected     deficit       sd
## 1 364132.31 -51814.3141 603.4338
## 2  46098.94  -2332.9402 214.7066
## 3  70484.75  -7466.7499 265.4896
## 4  28929.04  -1260.0429 170.0854
## 5  17908.48   -702.4761 133.8226

# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_9x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_9 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-other",
       title    = "Cumulative Deficit Hospitilization for other Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_9

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_9_adj <- combined_plot_9 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_9x_adj <- deficit_plot_9x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

cumulative_deficit_plot_9_adj <- cumulative_deficit_plot_9 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_9 <- combined_plot_9_adj /
                       deficit_plot_9x_adj /
                       cumulative_deficit_plot_9_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_9

ggsave("Other_combined_figure.png",
       plot = final_combined_plot_9,
       width = 11,       
       height = 10.5,    
       dpi = 300,
       bg = "white")

Different Approuch

library(ggplot2)
library(patchwork)
library(cowplot)
library(scales)

# Common scale settings
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(x_limits[1], x_limits[2], by = "1 year")
category_labels <- c("Cardiovascular", "Digestive", "Infectious", "Infectious Respiratory",
                     "Injury", "Pneumonia", "Renal", "Respiratory", "Other")

# ----- CLEANERS -----
clean_combined <- function(plots, labels) {
  mapply(function(p, label) {
    p +
      labs(title = label, x = NULL) +
      scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y") +
      theme_minimal() +
      theme(legend.position = "none",
            plot.title = element_text(hjust = 0.5, size = 11),
            axis.text.x = element_text(size = 8))
  }, plots, labels, SIMPLIFY = FALSE)
}

clean_deficit <- function(plot_list, labels) {
  mapply(function(p, label) {
    p +
      labs(title = label, x = NULL) +
      scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y") +
      theme_minimal() +
      theme(plot.title = element_text(hjust = 0.5, size = 11),
            axis.text.x = element_text(size = 8))
  }, plot_list, labels, SIMPLIFY = FALSE)
}

clean_cumulative <- function(plots, labels) {
  mapply(function(p, label) {
    p +
      labs(title = label, x = NULL, y = NULL, subtitle = NULL) +
      scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y") +
      scale_y_continuous(labels = comma) +
      theme_minimal() +
      theme(plot.title = element_text(hjust = 0.5, size = 11),
            axis.text = element_text(size = 8))
  }, plots, labels, SIMPLIFY = FALSE)
}

# ----- APPLY CLEANING -----
combined_clean <- clean_combined(list(
  combined_plot_1, combined_plot_2, combined_plot_3,
  combined_plot_4, combined_plot_5, combined_plot_6,
  combined_plot_7, combined_plot_8, combined_plot_9
), category_labels)

deficit_clean <- clean_deficit(list(
  deficit_plot_x, deficit_plot_2x, deficit_plot_3x,
  deficit_plot_4x, deficit_plot_5x, deficit_plot_6x,
  deficit_plot_7x, deficit_plot_8x, deficit_plot_9x
), category_labels)

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

cumulative_clean <- clean_cumulative(list(
  cumulative_deficit_plot_1, cumulative_deficit_plot_2, cumulative_deficit_plot_3,
  cumulative_deficit_plot_4, cumulative_deficit_plot_5, cumulative_deficit_plot_6,
  cumulative_deficit_plot_7, cumulative_deficit_plot_8, cumulative_deficit_plot_9
), category_labels)

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

# ----- LEGEND FOR COMBINED -----
legend_shared <- get_legend(
  combined_plot_1 +
    theme(legend.position = "bottom",
          legend.title = element_text(size = 10),
          legend.text = element_text(size = 9))
)

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

## Warning in get_plot_component(plot, "guide-box"): Multiple components found;
## returning the first one. To return all, use `return_all = TRUE`.

# ----- FINAL COMBINED PLOT (balanced, readable) -----
combined_final <- wrap_plots(combined_clean, ncol = 3) +
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom")

combined_final <- combined_final +
  plot_annotation(
    title = "Weekly Observed vs Expected Hospitalizations by Cause",
    theme = theme(plot.title = element_text(hjust = 0.5, size = 15))
  )

ggsave("combined_all_categories.png",
       plot = combined_final,
       width = 12, height = 10, dpi = 300, bg = "white")

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

deficit_final <- wrap_plots(deficit_clean, ncol = 3) +
  plot_annotation(
    title = "Deviations from Expected Weekly Admissions",
    theme = theme(plot.title = element_text(hjust = 0.5, size = 14))
  )

ggsave("deficit_all_categories.png",
       plot = deficit_final,
       width = 11, height = 9.5, dpi = 300, bg = "white")


cumulative_grid <- wrap_plots(cumulative_clean, ncol = 3)

y_label <- ggdraw() +
  draw_label("Cumulative Deficit in Hospitalizations", angle = 90, size = 12)

cumulative_final <- plot_grid(
  y_label, cumulative_grid,
  rel_widths = c(0.04, 0.96),
  nrow = 1
)

cumulative_final <- cumulative_final +
  plot_annotation(
    title = "Cumulative Deficit in Hospital Admissions by Cause (2020–2024)",
    theme = theme(plot.title = element_text(hjust = 0.5, size = 15))
  )

ggsave("cumulative_all_categories.png",
       plot = cumulative_final,
       width = 12, height = 10, dpi = 300, bg = "white")

combined_final

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

deficit_final

cumulative_final

ggsave("deficit_plot_x.png", plot = deficit_plot_x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_2x.png", plot = deficit_plot_2x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_3x.png", plot = deficit_plot_3x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_4x.png", plot = deficit_plot_4x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_5x.png", plot = deficit_plot_5x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_6x.png", plot = deficit_plot_6x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_7x.png", plot = deficit_plot_7x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_8x.png", plot = deficit_plot_8x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_9x.png", plot = deficit_plot_9x, width = 8, height = 6, dpi = 300)

Cumulative deficit interactive plots

cumulative_deficit <- deficit_cumulative(fit_x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for Cardio Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-Cardio", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_2x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_2 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-Digestive",
       title    = "Cumulative Deficit Hospitilization for Digestive Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_2

cumulative_deficit <- deficit_cumulative(fit_2x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for Digestive Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-Digestive", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#3

cumulative_deficit <- deficit_cumulative(fit_3x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for Infectious Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-Infectious", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#4

cumulative_deficit <- deficit_cumulative(fit_4x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for infectious respiratory Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-infectious respiratory", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#5

cumulative_deficit <- deficit_cumulative(fit_5x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for injury Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-injury", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#6

cumulative_deficit <- deficit_cumulative(fit_6x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for pneumonia Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-pneumonia", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#7

cumulative_deficit <- deficit_cumulative(fit_7x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for renal Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-renal", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#8

cumulative_deficit <- deficit_cumulative(fit_8x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for respiratory Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-respiratory", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#9

cumulative_deficit <- deficit_cumulative(fit_9x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for other Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-other", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

ggsave("cumulative_deficit_plot_1.png", plot = cumulative_deficit_plot_1, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_2.png", plot = cumulative_deficit_plot_2, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_3.png", plot = cumulative_deficit_plot_3, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_4.png", plot = cumulative_deficit_plot_4, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_5.png", plot = cumulative_deficit_plot_5, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_6.png", plot = cumulative_deficit_plot_6, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_7.png", plot = cumulative_deficit_plot_7, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_8.png", plot = cumulative_deficit_plot_8, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_9.png", plot = cumulative_deficit_plot_9, width = 8, height = 6, dpi = 300)

ggsave("cumulative_deficit_plot_1.png", plot = cumulative_deficit_plot_1, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_2.png", plot = cumulative_deficit_plot_2, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_3.png", plot = cumulative_deficit_plot_3, width = 5, height = 3, dpi = 300)

ggsave("cumulative_deficit_plot_4.png", plot = cumulative_deficit_plot_4, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_5.png", plot = cumulative_deficit_plot_5, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_6.png", plot = cumulative_deficit_plot_6, width = 5, height = 3, dpi = 300)

ggsave("cumulative_deficit_plot_7.png", plot = cumulative_deficit_plot_7, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_8.png", plot = cumulative_deficit_plot_8, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_9.png", plot = cumulative_deficit_plot_9, width = 5, height = 3, dpi = 300)

# new 'category' column added
cardiovascular <- fit_x$detected_intervals %>%
  mutate(category = "cardiovascular")
digestive <- fit_2x$detected_intervals %>%
  mutate(category = "digestive")
infectious <- fit_3x$detected_intervals %>%
  mutate(category = "infectious")
infectious_respiratory <- fit_4x$detected_intervals %>%
  mutate(category = "infectious respiratory")
injury <- fit_5x$detected_intervals %>%
  mutate(category = "injury")
pneumonia <- fit_6x$detected_intervals %>%
  mutate(category = "pneumonia")
renal <- fit_7x$detected_intervals %>%
  mutate(category = "renal")
respiratory <- fit_8x$detected_intervals %>%
  mutate(category = "respiratory")
other <- fit_9x$detected_intervals %>%
  mutate(category = "other")

combined_intervals <- bind_rows(
  cardiovascular,
  digestive,
  infectious,
  infectious_respiratory,
  injury,
  pneumonia,
  renal,
  respiratory,
  other
)

combined_intervals <- combined_intervals %>%
  select(category, everything())

combined_intervals <- combined_intervals %>%
  mutate(across(c(where(is.numeric)), ~round(., 1)))

kable_table <- combined_intervals %>%
  kable(
    format = "html",
    caption = "Detected Intervals by Disease Category",
    align = c("l", rep("c", 9)),
    col.names = c(
      "Category", "Start Date", "End Date", "Observed Count Rate", 
      "Expected Count Rate", "SD Count Rate", "Observed Counts", 
      "Expected Counts", "Deficit", "SD"
    ),
    row.names = FALSE
  )

kable_table

Detected Intervals by Disease Category

Category	Start Date	End Date	Observed Count Rate	Expected Count Rate	SD Count Rate	Observed Counts	Expected Counts	Deficit	SD
cardiovascular	2019-09-23	2019-12-02	8.1	8.7	0.1	12062	12840.1	-778.1	113.3
cardiovascular	2020-01-20	2020-02-03	8.3	8.7	0.1	3343	3515.5	-172.5	59.3
cardiovascular	2020-03-02	2021-02-15	7.3	8.9	0.0	49893	60938.8	-11045.8	246.9
cardiovascular	2021-12-27	2022-01-31	8.3	9.0	0.1	6709	7290.1	-581.1	85.4
cardiovascular	2023-05-01	2023-05-01	9.7	9.6	0.3	1303	1289.0	14.0	35.9
digestive	2019-09-15	2019-10-13	4.0	4.3	0.1	2684	2869.1	-185.1	53.6
digestive	2020-03-01	2020-09-06	3.5	4.4	0.0	13267	16748.9	-3481.9	129.4
digestive	2020-10-11	2021-03-14	3.7	4.2	0.0	11504	12987.3	-1483.3	114.0
digestive	2021-12-05	2022-01-30	3.8	4.2	0.1	4556	5073.5	-517.5	71.2
digestive	2022-05-01	2022-05-22	4.4	4.6	0.1	2353	2502.5	-149.5	50.0
digestive	2022-09-18	2022-10-02	4.3	4.5	0.1	1722	1835.4	-113.4	42.8
infectious	2020-05-03	2020-08-02	2.4	2.7	0.0	4495	5124.4	-629.4	71.6
infectious	2020-10-04	2021-05-30	2.2	2.7	0.0	10517	12848.3	-2331.3	113.4
infectious	2022-01-23	2022-02-27	2.5	2.8	0.1	2044	2250.7	-206.7	47.4
infectious	2022-12-25	2023-01-22	2.6	2.8	0.1	1727	1898.7	-171.7	43.6
infectious	2023-06-11	2023-07-02	2.5	2.8	0.1	1357	1490.1	-133.1	38.6
infectious	2024-10-20	2024-12-15	2.5	2.8	0.0	3060	3354.2	-294.2	57.9
infectious respiratory	2019-10-13	2019-12-01	0.6	1.3	0.0	683	1355.1	-672.1	36.8
infectious respiratory	2020-03-29	2021-05-30	0.1	0.9	0.0	1036	7153.0	-6117.0	84.6
infectious respiratory	2021-10-31	2022-03-20	0.5	1.3	0.0	1333	3752.6	-2419.6	61.3
infectious respiratory	2023-01-15	2023-02-26	0.7	1.3	0.0	616	1184.6	-568.6	34.4
infectious respiratory	2023-04-09	2023-04-23	0.4	0.7	0.0	170	283.3	-113.3	16.8
infectious respiratory	2024-11-03	2024-11-10	0.6	1.0	0.1	173	270.3	-97.3	16.4
injury	2019-10-20	2019-11-24	6.8	7.3	0.1	5458	5880.9	-422.9	76.7
injury	2020-01-12	2020-09-13	6.3	7.6	0.0	30779	36951.6	-6172.6	192.2
injury	2020-10-04	2021-02-07	6.6	7.5	0.1	16830	19168.3	-2338.3	138.4
injury	2021-03-07	2021-04-11	7.1	7.6	0.1	5752	6123.8	-371.8	78.3
injury	2021-05-02	2021-05-16	7.6	7.9	0.1	3070	3208.8	-138.8	56.6
injury	2021-07-25	2021-08-08	8.0	8.4	0.1	3223	3381.7	-158.7	58.2
injury	2021-12-19	2022-01-16	7.1	7.7	0.1	4807	5171.7	-364.7	71.9
injury	2022-07-31	2022-08-21	8.3	8.7	0.1	4453	4668.8	-215.8	68.3
pneumonia	2019-10-27	2019-12-01	2.5	3.0	0.1	2046	2397.5	-351.5	49.0
pneumonia	2020-04-26	2021-07-04	1.5	2.6	0.0	12636	22211.3	-9575.3	149.0
pneumonia	2021-09-05	2022-05-22	1.9	2.8	0.0	9964	14396.8	-4432.8	120.0
pneumonia	2022-07-17	2022-08-28	1.7	2.0	0.0	1597	1931.8	-334.8	44.0
pneumonia	2022-11-06	2022-11-13	2.5	2.9	0.1	683	786.8	-103.8	28.1
pneumonia	2023-01-22	2023-02-19	2.4	3.0	0.1	1621	2005.2	-384.2	44.8
pneumonia	2023-04-23	2023-05-21	2.3	2.7	0.1	1545	1792.5	-247.5	42.3
pneumonia	2024-01-28	2024-02-11	2.6	3.0	0.1	1058	1201.7	-143.7	34.7
renal	2019-11-10	2019-12-01	1.9	2.1	0.1	1031	1114.3	-83.3	33.4
renal	2020-03-01	2020-07-12	1.8	2.3	0.0	4844	6219.3	-1375.3	78.9
renal	2020-09-06	2020-10-04	2.2	2.3	0.1	1461	1566.4	-105.4	39.6
renal	2020-11-01	2021-02-28	2.0	2.2	0.0	4765	5376.3	-611.3	73.3
renal	2021-12-26	2022-01-02	2.2	2.3	0.1	596	627.6	-31.6	25.1
renal	2022-05-15	2022-05-29	2.4	2.6	0.1	955	1053.3	-98.3	32.5
respiratory	2019-07-21	2019-07-21	3.0	3.4	0.2	401	453.3	-52.3	21.3
respiratory	2019-10-27	2019-11-24	3.9	4.3	0.1	2639	2911.5	-272.5	54.0
respiratory	2020-03-08	2021-05-30	2.7	4.2	0.0	23417	36478.4	-13061.4	191.0
respiratory	2021-11-21	2022-03-27	3.6	4.4	0.0	9293	11363.2	-2070.2	106.6
respiratory	2023-04-23	2023-05-07	4.2	4.7	0.1	1691	1897.1	-206.1	43.6
other	2020-03-01	2021-05-23	35.7	41.6	0.1	312318	364132.3	-51814.3	603.4
other	2021-08-22	2021-10-10	40.6	42.8	0.2	43766	46098.9	-2332.9	214.7
other	2021-11-21	2022-02-13	36.0	40.3	0.2	63018	70484.7	-7466.7	265.5
other	2022-03-06	2022-04-03	41.1	43.0	0.3	27669	28929.0	-1260.0	170.1
other	2022-04-24	2022-05-08	42.6	44.3	0.3	17206	17908.5	-702.5	133.8

library(kableExtra)
html_table <- combined_intervals %>%
  kable("html",
        escape = FALSE,
        caption = "Detected Intervals by Disease Category",
        align = c("l", rep("c", ncol(combined_intervals) - 1)),
        col.names = c(
          "Category", "Start Date", "End Date", "Observed Count Rate", 
          "Expected Count Rate", "SD Count Rate", "Observed Counts", 
          "Expected Counts", "Deficit", "SD"
        )) %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                position = "center")

html_table

Detected Intervals by Disease Category

	Category	Start Date	End Date	Observed Count Rate	Expected Count Rate	SD Count Rate	Observed Counts	Expected Counts	Deficit	SD
1…1	cardiovascular	2019-09-23	2019-12-02	8.1	8.7	0.1	12062	12840.1	-778.1	113.3
2…2	cardiovascular	2020-01-20	2020-02-03	8.3	8.7	0.1	3343	3515.5	-172.5	59.3
3…3	cardiovascular	2020-03-02	2021-02-15	7.3	8.9	0.0	49893	60938.8	-11045.8	246.9
4…4	cardiovascular	2021-12-27	2022-01-31	8.3	9.0	0.1	6709	7290.1	-581.1	85.4
5…5	cardiovascular	2023-05-01	2023-05-01	9.7	9.6	0.3	1303	1289.0	14.0	35.9
1…6	digestive	2019-09-15	2019-10-13	4.0	4.3	0.1	2684	2869.1	-185.1	53.6
2…7	digestive	2020-03-01	2020-09-06	3.5	4.4	0.0	13267	16748.9	-3481.9	129.4
3…8	digestive	2020-10-11	2021-03-14	3.7	4.2	0.0	11504	12987.3	-1483.3	114.0
4…9	digestive	2021-12-05	2022-01-30	3.8	4.2	0.1	4556	5073.5	-517.5	71.2
5…10	digestive	2022-05-01	2022-05-22	4.4	4.6	0.1	2353	2502.5	-149.5	50.0
6…11	digestive	2022-09-18	2022-10-02	4.3	4.5	0.1	1722	1835.4	-113.4	42.8
1…12	infectious	2020-05-03	2020-08-02	2.4	2.7	0.0	4495	5124.4	-629.4	71.6
2…13	infectious	2020-10-04	2021-05-30	2.2	2.7	0.0	10517	12848.3	-2331.3	113.4
3…14	infectious	2022-01-23	2022-02-27	2.5	2.8	0.1	2044	2250.7	-206.7	47.4
4…15	infectious	2022-12-25	2023-01-22	2.6	2.8	0.1	1727	1898.7	-171.7	43.6
5…16	infectious	2023-06-11	2023-07-02	2.5	2.8	0.1	1357	1490.1	-133.1	38.6
6…17	infectious	2024-10-20	2024-12-15	2.5	2.8	0.0	3060	3354.2	-294.2	57.9
1…18	infectious respiratory	2019-10-13	2019-12-01	0.6	1.3	0.0	683	1355.1	-672.1	36.8
2…19	infectious respiratory	2020-03-29	2021-05-30	0.1	0.9	0.0	1036	7153.0	-6117.0	84.6
3…20	infectious respiratory	2021-10-31	2022-03-20	0.5	1.3	0.0	1333	3752.6	-2419.6	61.3
4…21	infectious respiratory	2023-01-15	2023-02-26	0.7	1.3	0.0	616	1184.6	-568.6	34.4
5…22	infectious respiratory	2023-04-09	2023-04-23	0.4	0.7	0.0	170	283.3	-113.3	16.8
6…23	infectious respiratory	2024-11-03	2024-11-10	0.6	1.0	0.1	173	270.3	-97.3	16.4
1…24	injury	2019-10-20	2019-11-24	6.8	7.3	0.1	5458	5880.9	-422.9	76.7
2…25	injury	2020-01-12	2020-09-13	6.3	7.6	0.0	30779	36951.6	-6172.6	192.2
3…26	injury	2020-10-04	2021-02-07	6.6	7.5	0.1	16830	19168.3	-2338.3	138.4
4…27	injury	2021-03-07	2021-04-11	7.1	7.6	0.1	5752	6123.8	-371.8	78.3
5…28	injury	2021-05-02	2021-05-16	7.6	7.9	0.1	3070	3208.8	-138.8	56.6
6…29	injury	2021-07-25	2021-08-08	8.0	8.4	0.1	3223	3381.7	-158.7	58.2
7…30	injury	2021-12-19	2022-01-16	7.1	7.7	0.1	4807	5171.7	-364.7	71.9
8…31	injury	2022-07-31	2022-08-21	8.3	8.7	0.1	4453	4668.8	-215.8	68.3
1…32	pneumonia	2019-10-27	2019-12-01	2.5	3.0	0.1	2046	2397.5	-351.5	49.0
2…33	pneumonia	2020-04-26	2021-07-04	1.5	2.6	0.0	12636	22211.3	-9575.3	149.0
3…34	pneumonia	2021-09-05	2022-05-22	1.9	2.8	0.0	9964	14396.8	-4432.8	120.0
4…35	pneumonia	2022-07-17	2022-08-28	1.7	2.0	0.0	1597	1931.8	-334.8	44.0
5…36	pneumonia	2022-11-06	2022-11-13	2.5	2.9	0.1	683	786.8	-103.8	28.1
6…37	pneumonia	2023-01-22	2023-02-19	2.4	3.0	0.1	1621	2005.2	-384.2	44.8
7…38	pneumonia	2023-04-23	2023-05-21	2.3	2.7	0.1	1545	1792.5	-247.5	42.3
8…39	pneumonia	2024-01-28	2024-02-11	2.6	3.0	0.1	1058	1201.7	-143.7	34.7
1…40	renal	2019-11-10	2019-12-01	1.9	2.1	0.1	1031	1114.3	-83.3	33.4
2…41	renal	2020-03-01	2020-07-12	1.8	2.3	0.0	4844	6219.3	-1375.3	78.9
3…42	renal	2020-09-06	2020-10-04	2.2	2.3	0.1	1461	1566.4	-105.4	39.6
4…43	renal	2020-11-01	2021-02-28	2.0	2.2	0.0	4765	5376.3	-611.3	73.3
5…44	renal	2021-12-26	2022-01-02	2.2	2.3	0.1	596	627.6	-31.6	25.1
6…45	renal	2022-05-15	2022-05-29	2.4	2.6	0.1	955	1053.3	-98.3	32.5
1…46	respiratory	2019-07-21	2019-07-21	3.0	3.4	0.2	401	453.3	-52.3	21.3
2…47	respiratory	2019-10-27	2019-11-24	3.9	4.3	0.1	2639	2911.5	-272.5	54.0
3…48	respiratory	2020-03-08	2021-05-30	2.7	4.2	0.0	23417	36478.4	-13061.4	191.0
4…49	respiratory	2021-11-21	2022-03-27	3.6	4.4	0.0	9293	11363.2	-2070.2	106.6
5…50	respiratory	2023-04-23	2023-05-07	4.2	4.7	0.1	1691	1897.1	-206.1	43.6
1…51	other	2020-03-01	2021-05-23	35.7	41.6	0.1	312318	364132.3	-51814.3	603.4
2…52	other	2021-08-22	2021-10-10	40.6	42.8	0.2	43766	46098.9	-2332.9	214.7
3…53	other	2021-11-21	2022-02-13	36.0	40.3	0.2	63018	70484.7	-7466.7	265.5
4…54	other	2022-03-06	2022-04-03	41.1	43.0	0.3	27669	28929.0	-1260.0	170.1
5…55	other	2022-04-24	2022-05-08	42.6	44.3	0.3	17206	17908.5	-702.5	133.8

save_kable(html_table, "combined_intervals_table.html")
webshot("combined_intervals_table.html", "combined_intervals_table.png", vwidth = 1200, vheight = 800)

combined_intervals <- combined_intervals %>%
  mutate(across(where(is.numeric), ~ round(., 1)))

table_1 <- combined_intervals %>%
  slice(1:5) %>%
  kbl(
    format = "html",
    caption = "Detected Intervals by Category (1 of 2)",
    align = c("l", rep("c", 9)),
    col.names = c(
      "Category", "Start Date", "End Date", "Observed Count Rate", 
      "Expected Count Rate", "SD Count Rate", "Observed Counts", 
      "Expected Counts", "Deficit", "SD"
    )
  ) %>%
  kable_styling(full_width = FALSE)
save_kable(table_1, file = "detected_intervals_1.png", zoom = 3)

## save_kable will have the best result with magick installed.

table_2 <- combined_intervals %>%
  slice(6:9) %>%
  kbl(
    format = "html",
    caption = "Detected Intervals by Category (2 of 2)",
    align = c("l", rep("c", 9)),
    col.names = c(
      "Category", "Start Date", "End Date", "Observed Count Rate", 
      "Expected Count Rate", "SD Count Rate", "Observed Counts", 
      "Expected Counts", "Deficit", "SD"
    )
  ) %>%
  kable_styling(full_width = FALSE)
save_kable(table_2, file = "detected_intervals_2.png", zoom = 3)

## save_kable will have the best result with magick installed.

table_2 

Detected Intervals by Category (2 of 2)

Category	Start Date	End Date	Observed Count Rate	Expected Count Rate	SD Count Rate	Observed Counts	Expected Counts	Deficit	SD
digestive	2019-09-15	2019-10-13	4.0	4.3	0.1	2684	2869.1	-185.1	53.6
digestive	2020-03-01	2020-09-06	3.5	4.4	0.0	13267	16748.9	-3481.9	129.4
digestive	2020-10-11	2021-03-14	3.7	4.2	0.0	11504	12987.3	-1483.3	114.0
digestive	2021-12-05	2022-01-30	3.8	4.2	0.1	4556	5073.5	-517.5	71.2

Further Analysis

library(ggplot2)
library(scales)

total_deficits <- data.frame(
  category = c("Cardiovascular", "Digestive", "Infectious", "Infectious Respiratory",
               "Injury", "Pneumonia", "Renal", "Respiratory", "Other"),
  total_deficit = c(
    sum(cumulative_deficit_plot_1$data$observed, na.rm = TRUE),
    sum(cumulative_deficit_plot_2$data$observed, na.rm = TRUE),
    sum(cumulative_deficit_plot_3$data$observed, na.rm = TRUE),
    sum(cumulative_deficit_plot_4$data$observed, na.rm = TRUE),
    sum(cumulative_deficit_plot_5$data$observed, na.rm = TRUE),
    sum(cumulative_deficit_plot_6$data$observed, na.rm = TRUE),
    sum(cumulative_deficit_plot_7$data$observed, na.rm = TRUE),
    sum(cumulative_deficit_plot_8$data$observed, na.rm = TRUE),
    sum(cumulative_deficit_plot_9$data$observed, na.rm = TRUE)
  )
)

bar_plot <- ggplot(total_deficits, aes(x = reorder(category, total_deficit), y = total_deficit)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  scale_y_continuous(labels = comma) +
  labs(title = "Total Cumulative Deficits by Category (2020–2024)",
       x = "ICD-10 Category",
       y = "Total Cumulative Deficit in Hospitalizations") +
  coord_flip() +
  theme_minimal(base_size = 12)
bar_plot

ggsave("Total_Cumulative_Deficit_by_Category.png",
       plot = bar_plot,
       width = 8.5,
       height = 6,
       dpi = 300,
       bg = "white")

cardio_deficit     <- tail(deficit_cumulative(fit_x,   start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
digestive_deficit  <- tail(deficit_cumulative(fit_2x,  start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
infectious_deficit <- tail(deficit_cumulative(fit_3x,  start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
inf_resp_deficit   <- tail(deficit_cumulative(fit_4x,  start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
injury_deficit     <- tail(deficit_cumulative(fit_5x,  start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
pneumonia_deficit  <- tail(deficit_cumulative(fit_6x,  start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
renal_deficit      <- tail(deficit_cumulative(fit_7x,  start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
respiratory_deficit<- tail(deficit_cumulative(fit_8x,  start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
other_deficit      <- tail(deficit_cumulative(fit_9x,  start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed

deficit_summary <- tibble(
  Category = c("Cardiovascular", "Digestive", "Infectious", "Infectious Respiratory", 
               "Injury", "Pneumonia", "Renal", "Respiratory", "Other"),
  Cumulative_Deficit = c(cardio_deficit, digestive_deficit, infectious_deficit,
                         inf_resp_deficit, injury_deficit, pneumonia_deficit,
                         renal_deficit, respiratory_deficit, other_deficit)
)

bar_plot2 <- ggplot(deficit_summary, aes(x = reorder(Category, -Cumulative_Deficit), 
                            y = Cumulative_Deficit)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  geom_text(aes(label = scales::comma(round(Cumulative_Deficit))), 
            vjust = -0.5, size = 3) +
  scale_y_continuous(labels = scales::comma) +
  labs(
    title = "Cumulative Deficit in Hospitalizations by Category",
    subtitle = "Massachusetts, March 2020 – December 2024",
    x = "ICD-10 Category",
    y = "Cumulative Deficit (Observed - Expected)"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
bar_plot2

ggsave("Cumulative Deficit in Hospitalizations by Category.png",
       plot = bar_plot2,
       width = 8.5,
       height = 6,
       dpi = 300,
       bg = "white")

ICD-10 codes table

library(tidyverse)
library(readxl)
library(gt)
library(stringr)

# Load and clean data
icd10categories <- read_excel("icd10categories_updated_321.xlsx") %>%
  mutate(icd10code = toupper(trimws(icd10code)))

# Function to generate code ranges and extract sort keys
get_code_range_with_sort <- function(codes) {
  codes <- sort(unique(codes))
  letter <- substr(codes, 1, 1)
  num <- as.integer(str_extract(substr(codes, 2, 4), "\\d+"))
  df <- tibble(code = codes, letter = letter, num = num)
  
  ranges <- df %>%
    group_by(letter) %>%
    summarize(
      min_num = min(num),
      max_num = max(num),
      .groups = "drop"
    ) %>%
    mutate(
      range = ifelse(min_num == max_num,
                     paste0(letter, sprintf("%02d", min_num)),
                     paste0(letter, sprintf("%02d", min_num), "–", letter, sprintf("%02d", max_num))),
      sort_key = min_num + (utf8ToInt(letter) * 1000) # sort letter first, then number
    )
  
  list(
    range = paste(ranges$range, collapse = ", "),
    sort_key = min(ranges$sort_key)
  )
}

# Apply function and unpack results
icd10_summary <- icd10categories %>%
  group_by(category_0321update) %>%
  summarize(
    `Number of ICD-10 codes` = n_distinct(icd10code),
    .groups = "drop",
    temp = list(get_code_range_with_sort(icd10code))
  ) %>%
  mutate(
    `ICD-10 code groups` = map_chr(temp, "range"),
    sort_key = map_dbl(temp, "sort_key")
  ) %>%
  select(-temp) %>%
  rename(category = category_0321update) %>%
  arrange(sort_key)  # Sort by ICD-10 code order

## Warning: There were 3 warnings in `summarize()`.
## The first warning was:
## ℹ In argument: `sort_key = min_num + (utf8ToInt(letter) * 1000)`.
## ℹ In group 0: .
## Caused by warning:
## ! There was 1 warning in `mutate()`.
## ℹ In argument: `sort_key = min_num + (utf8ToInt(letter) * 1000)`.
## Caused by warning in `utf8ToInt()`:
## ! argument should be a character vector of length 1
## all but the first element will be ignored
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 2 remaining warnings.

# Create GT table
gt_table3 <- icd10_summary %>%
  select(-sort_key) %>%
  gt() %>%
  tab_header(
    title = md("**Table. ICD-10 Code Groupings and Counts by Category**")
  ) %>%
  cols_width(
    `ICD-10 code groups` ~ px(400)
  ) %>%
  tab_options(table.width = pct(100))
gt_table3

Table. ICD-10 Code Groupings and Counts by Category
category	Number of ICD-10 codes	ICD-10 code groups
infectious	199	A00–A99, B01–B99
cardiovascular	100	I00–I99
infectious respiratory	13	J00–J22
pneumonia	7	J12–J18
respiratory	78	J19–J99
digestive	25	K20–K92
injury	342	M48–M97, R45, S00–S99, T00–T88, V00–V98, W00–W88, X00–X99, Y00–Y99
renal	29	N00–N28
other	75	O00–O99, Z31–Z33
covid	1	U07

# Save the GT table
gtsave(gt_table3, filename = "icd10_category_groupings_continuous.png")

library(tidyverse)
library(readxl)
library(gt)
library(stringr)

# Load and clean data
icd10categories <- read_excel("icd10categories_updated_321.xlsx") %>%
  mutate(icd10code = toupper(trimws(icd10code)))

# Modified condensation function with relaxed grouping logic
condense_icd10_codes <- function(codes, max_gap = 2) {
  codes <- sort(unique(codes))
  letter_part <- substr(codes, 1, 1)
  num_part <- as.integer(str_extract(substr(codes, 2, 4), "\\d+"))
  
  df <- tibble(code = codes, letter = letter_part, num = num_part) %>%
    arrange(letter, num)
  
  # Relax grouping: allow small gaps (max_gap)
  df <- df %>%
    arrange(letter, num) %>%
    group_by(letter) %>%
    mutate(group = cumsum(c(TRUE, diff(num) > max_gap))) %>%
    group_by(letter, group) %>%
    summarize(
      range = if (n() == 1) first(code) else paste0(first(code), "–", last(code)),
      sort_val = first(num) + utf8ToInt(first(letter)) * 1000,
      .groups = "drop"
    )
  
  list(
    range = paste(df$range, collapse = ", "),
    sort_key = min(df$sort_val)
  )
}

# Summarize and apply function
icd10_summary <- icd10categories %>%
  group_by(category_0321update) %>%
  summarize(
    `Number of ICD-10 codes` = n_distinct(icd10code),
    temp = list(condense_icd10_codes(icd10code, max_gap = 2)),
    .groups = "drop"
  ) %>%
  mutate(
    `ICD-10 code groups` = map_chr(temp, "range"),
    sort_key = map_dbl(temp, "sort_key")
  ) %>%
  rename(category = category_0321update) %>%
  arrange(sort_key)  # Alphabetical ICD-10 ordering

# Create GT table
gt_table <- icd10_summary %>%
  select(-sort_key, -temp) %>%
  gt() %>%
  tab_header(
    title = md("**Table 1. ICD-10 Code Groupings and Counts by Category**")
  ) %>%
  cols_width(
    `ICD-10 code groups` ~ px(500)
  ) %>%
  tab_options(table.width = pct(100))
gt_table

Table 1. ICD-10 Code Groupings and Counts by Category
category	Number of ICD-10 codes	ICD-10 code groups
infectious	199	A00–A99, B01–B99
cardiovascular	100	I00–I99
infectious respiratory	13	J00–J06, J09–J11, J20–J22
pneumonia	7	J12–J18
respiratory	78	J19, J23–J99
digestive	25	K20–K31, K50–K52, K55–K63, K92
injury	342	M48, M80, M84, M96–M97, R45, S00–S99, T00–T88, V00–V05, V09–V15, V18–V29, V40–V49, V53–V74, V79–V83, V86–V98, W00–W01, W04–W22, W25–W31, W34, W39–W40, W44–W46, W49–W57, W61, W67–W69, W74, W86–W88, X00–X02, X08, X14, X19, X30, X58, X74, X78–X80, X83, X95, X99, Y00, Y03–Y04, Y07–Y09, Y21, Y24, Y27–Y32, Y63, Y69–Y73, Y79, Y82–Y84, Y90–Y95, Y99
renal	29	N00–N28
other	75	O00–O04, O07–O16, O20–O36, O40–O48, O60–O77, O80–O82, O85–O94, O98–O99, Z31–Z33
covid	1	U07

# Save table
gtsave(gt_table, filename = "icd10_category_groupings.png")